CERN/TH-99-376 
CPT-99/PE.3916 

LAPTH-Conf-771 /99 



O 
O 

^ : A numerical treatment of Neuberger's lattice 

; Dirac operator* 

Pilar Hernandez| Karl Jansen^ 

> 
00 

O ■ CERN, 1211 Geneva 23, Switzerland 

O 
O 

^ ! and Laurent Lellouch^ 

> 

g^' LAPTH, Chemin de Bellevue, B.P. 110, 

^ ■ 

F74941 Annecy-le-Vieux Cedex, France 



X 



February 1, 2008 

Abstract 

We describe in some detail our numerical treatment of Neuberger's lat- 
tice Dirac operator as implemented in a practical application. We discuss 
the improvements we have found to accelerate the numerical computations 
and give an estimate of the expense when using this operator in practice. 
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1 Lattice formulation of QCD 



Today, we believe that the world of quarks and gluons is described theoreti- 
cally by quantum chromodynamics (QCD). This model shows a number of non- 
perturbative aspects that cannot be adequately addressed by approximation schemes 
such as perturbation theory. The only way to evaluate QCD, addressing both its 
perturbative and non-perturbative aspects at the same time, is lattice QCD. In 
this approach the theory is put on a 4-dimensional Euclidean space-time lattice 
of finite physical length L, with a no n- vanishing value of the lattice spacing a. 
Having only a finite number of grid points, physical quantities can be computed 
numerically by solving a high- dimensional integral by Monte Carlo methods, mak- 
ing use of importance sampling. 

The introduction of a lattice spacing regularizes the theory and is an intermediate 
step in the computation of physical observables. Eventually, the regularization 
has to be removed and the value of the lattice spacing has to be sent to zero to 
reach the target theory, i.e. continuum QCD. In fact, in conventional formulations 
of lattice QCD the introduction of the lattice spacing renders the theory 
on the lattice somewhat different from the continuum analogue and a number 
of properties of the continuum theory are only very difficult and cumbersome 
to establish in the lattice regularized theory. One of the main reasons for this 
difficulty is that in conventional lattice QCD the regularization breaks a particular 
symmetry of the continuum theory, which plays a most important role there, 
namely chiral symmetry. 

However, the last few years have seen a major breakthrough in that we now have 
formulations of lattice QCD that have an exact lattice chiral symmetry [Q. In 
this approach, many properties of continuum QCD are preserved even at a non- 
vanishing value of the lattice spacing 0, ^ ||, ||, ^. This development followed 
the rediscovery [0 of the so-called Ginsparg-Wilson (GW) relation [Q which is 
fulfilled by any operator with the exact lattice chiral symmetry of It is not the 
aim of this contribution to discuss the physics consequences of the GW relation. 
We have to refer the interested reader to reviews |TU| about these topics. Here 
we would like to discuss the numerical treatment of a particular lattice operator 
that satisfies the GW relation, namely Neuberger's solution . This solution has 
a complicated structure and is challenging to implement numerically. Thus, the 
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large theoretical advantage of an operator satisfying the GW relation must be 
weighed against the very demanding computational effort required to implement 
it. 

This contribution is organized as follows. After discussing Neuberger's lattice 
Dirac operator we want to show how we evaluated the operator in our practi- 



cal application [|lT| and what kind of improvements we found to accelerate the 
numerical computations. For alternative ideas for improvements, see the contri- 
butions of H. Neuberger |T2| and A. Borici to this workshop. We finally give 
some estimates of the computational expense of using Neuberger's operator. 



2 Neuberger's lattice Dirac operator 

The operator we have used acts on fields (complex vectors) where x = 

{xq, xi, X2, X3) and the x^, /i = 0, 1, 2, 3, are integer numbers denoting a 4-dimensional 
grid point in a lattice of size N"^ with N = L/a. The fields carry in addition 
a "colour" index a = 1,2,3 as well as a "Dirac" index i = 1,2,3,4. Hence, $ is 
a A^^ ■ 3 ■ 4 complex vector. 

In order to reach the expression for Neuberger's operator we first introduce the 
matrix A 

A = l + s-^{^,{V; + W,)-aW;V,} , (1) 

where and V* are the nearest-neighbour forward and backward derivatives, 
the precise definition of which can be found in the Appendix. The parameter 
s is to be taken in the range of |s| < 1 and serves to optimize the localization 



properties [|1^ of Neuberger's operator, which is then given by 



D = -{l-A{A^Ay'^']. (2) 



a 

Through the appearance of the square root in eq. @, all points on the lattice are 
connected with one another, giving rise to a very complicated, multi-neighbour 
action. However, the application of D to a vector $ will only contain applications 
of A or A'^A on this vector. Since these matrices are sparse, as only nearest- 
neighbour interactions are involved, we will never have to store the whole matrix. 
In the computation of physical quantities, the inverse of D, applied to a given 
vector, is generically needed. Hence one faces the problem of having to compute 
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a vector X = D'^rj, with rj a prescribed vector (the "source") as required by the 
particular problem under investigation. Fortunately, a number of efficiently work- 
ing algorithms for computing X = D'^rj are known, such as conjugate gradient. 



BiCGstab, or variants thereof |T6[. In conventional approaches to lattice QCD an 
operator D is used that is very similar to the matrix A in eq. (|l]). Computing the 
vector X = D~^ri requires a number njter of iterations of some particular method, 
say BiCGstab. Employing Neuberger's operator D in computing X = D~^ri, it 
turns out that the number of iterations needed is of the same order of magnitude 
as when using D. At the same time, in each of these iterations, the square root 
has to be evaluated. When this is done by some polynomial approximation, it is 
found that the required degree of this polynomial is roughly of the same order as 
the number of iterations needed for computing the vector X. Hence, with respect 
to the conventional case, the numerical effort is squared and the price to pay for 
using the operator D is high. 

On the other hand, any solution of the Ginsparg-Wilson relation gives us a tool 
by which particular problems in lattice QCD can be studied, which would be 
extremely hard to address with conventional approaches. It is for these cases 
that the large numerical effort is justified, but clearly, we would like to have 
clever ideas coming from areas such as Applied Mathematics, to decrease the 
numerical expense or even overcome this bottleneck. 



3 Approximation of (^^^) 



-1/2 



For computing the square root that appears in eq. (|^), we have chosen a Cheby- 
shev approximation |T3| by constructing a polynomial Pn,e{x) of degree n, which 



has an exponential convergence rate in the interval x G [e, 1]. Outside this inter- 
val, convergence is still found but it will not be exponential. The advantages of 
using this Chebyshev approximation are the well-controlled exponential fit accu- 
racy as well as the possibility of having numerically very stable recursion relations 
?| to construct the polynomial, allowing for large degrees. In order to have an 



optimal approximation, it is desirable to know the lowest and the highest eigen- 
value of A^A. A typical example of the eigenvalues of A^A is shown in fig. 
where we show the 11 lowest eigenvalues as obtained on a number of configu- 
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Figure 1: Monte Carlo time evolution of the eleven lowest eigenvalues of A^A at 
f3 = 5.85. The lowest eigenvalue for each configuration is the open square. 



rations using the Ritz functional method There is a wide spread and very 
low-lying eigenvalues appear. Choosing e to be the value of the lowest of these 
eigenvalues would result in a huge degree n of the polynomial Pn^e- We therefore 
computed 0(10) lowest-lying eigenvalues of A'^A as well as their eigenfunctions 
and projected them out of the matrix A'^A. The approximation is then only per- 
formed for the matrix with a reduced condition number, resulting in a substantial 
decrease of the degree of the polynomial. In addition, we computed the highest 
eigenvalue of A^A and normalized the matrix A such that ||v4'''yl|| < 1. 



Since our work [|lT|, aiming at the physical question of spontaneous chiral sym- 
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metry breaking in lattice QCD, has been one of the first of its kind, we wanted 
to exclude possible systematic errors and demanded a very high precision for the 
approximation to the square root: 

||X - Pn,e{A^ A)A^ APnM^ A)X\\'^ /\\2X\\'^ < 10-^^ (3) 

where X is a gaussian random vector. In our practical applications we fixed this 
precision beforehand and set e to be the 11th lowest eigenvalue of A'^A. This then 
determines the degree of the polynomial n and hence our approximation D„ to 
the exact Neuberger operator D. We checked that the precision we required for 
the approximation of the square root is directly related to the precision by which 
the GW relation itself is fulfilled. Choosing n such that the accuracy in eq. (^) 
is reached results in 

II [75D„ + D„,75 - /^„75^n] X|| VllXf ^ 10-^^ . (4) 

In addition, we find that the deviations from the exact GW relation decrease 
exponentially fast with increasing n. 

4 The inverse of Neuberger's operator 

As mentioned above, in physics applications a vector D^^rj has to be computed, 
with 1] a prescribed source vector. Not only is the computation of this vector very 
costly, there also appears to be a conceptual problem: in inspecting the lowest 
eigenvalue of Dj^Dn, very small eigenvalues are often found as shown in fig. |^. 
These very small eigenvalues belong to a given chiral sector of the theory, i.e. 
their corresponding eigenfunctions x s^i's eigenfunctions of 75 with 75% = ±x- 
fact, these modes play an important physical role as they are associated with 
topological sectors of the theory 0, ||, ||, ^. 

As far as the practical applications are concerned, it is clear that in the presence 
of such a small eigenvalue, the inversion of Dn will be very costly, as the condition 
number of the problem is then very high. In order to address this problem, we 
followed two strategies: 

(i) We compute the lowest eigenvalue of Dj^Dn and its eigenfunction (using 
again the Ritz functional method [|I^) and if it is a zero mode -in which 
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Figure 2: Monte Carlo time evolution of the lowest eigenvalue of D^Dn- The 
eigenvalues belong to given chiral sectors of the theory denoted as ±x ior chirality 
plus (full squares) and minus (open squares). Data are obtained at (3 — 5.85 
choosing s — 0.6. Whenever there is a zero mode of Dj^Dn, the value of the 
lowest eigenvalue is set to 10~^. 
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case it is also a zero mode of we project this mode out of Dn and 
invert only the reduced matrix; this is then well conditioned, as the very 
small eigenvalues appear to be isolated. In this strategy, the knowledge of 
the eigenfunction must be very precise and an accuracy of approximating 
the square root as indicated in eq. (H) is mandatory. 

{ii) Again we determine the lowest eigenvalue of Dl^Dn and the chirality of the 
corresponding zero mode, if there is any. We then make use of the fact that 
Dl^Dn commutes with 75. This allows us to perform the inversion in the 
chiral sector without zero modes. In this strategy, the accuracy demanded in 



eq. (g) could be relaxed and this strategy, which essentially follows ref. |[T8[| , 
is in general much less expensive than following strategy (i). 

However, even adopting strategy [ii), solving the system DnX = rj is still costly. 
We therefore tried two ways of improving on this. We first note that instead of 
solving 



X = ri (5) 



we can equally well solve 



X = A^ri. (6) 



In practice, however, we found no real advantage in using the formulation of 
eq. (^). We have further considered two acceleration schemes. 

Scheme (a) 

We choose two different polynomials (now approximating V A^A and not the 
inverse) Pn,e and Pm,e, m < n, such that 

Pn,e = Pm,e + A (7) 

with A a "small" correction. Then we have 

"i + (At - p^;) a] (At - p^,,) . (8) 

This leads us to the following procedure of solving DnX = rj: 
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(1) first solve 



(At - P„,,) Y 



(9) 



(2) then solve 



[A^ -Pm,e)X^ = 7] + AY ■ (10) 

(3) use Xo as a starting vector to finally solve 



(At - P„,,) X = 7]. (11) 

The generation of the starting vector Xq in steps (1) and (2) is only a small 
overhead. In fig. ^ we plot the relative residuum eg^^p = ||-DnX — r7||^/||X||^ as a 
function of the number of applications of D^- In this case n = 100 and m = 30. 
We show the number of applications of the matrix Dn for the case of a random 
starting vector (dotted line) and the case where Xq was generated according to 
the above procedure (solid line). The gain is of approximately a factor of two. 

Scheme (b) 

In the second approach, we use a sequence of polynomials to solve -D„X = r]. To 
this end we first solve 

[A^-P^,,,]X, = A^r] (12) 

by choosing a polynomial Pmi,e and a stopping criterion for the solver e^^p such 
that 

mi < n, e£p > estop • (13) 

The value of e^^l^ is chosen such that it is roughly of the same order of magnitude 
as the error that the polynomial of degree rui itself induces. The solution Xi 
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Figure 3: The residuum as a function of the number of appHcations of the matrix 
Dji- The dotted hne corresponds to a random starting vector. The sohd hne to 
a starting vector generated following scheme (a). 
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is then used as a starting vector for the next equation, employing a polynomial 
Pm2,€ and stopping criterion el^^p with 

mi < ma < n, eJt^p > e^ip > Cstop ■ (14) 

This procedure is then repeated until we reach the desired polynomial Pn,t and 
stopping criterion estop to solve the real equation 



[A^ - P„,,] X = A^7j. (15) 
As for scheme (a), we gain a factor of about two in the numerical effort. We 



finally remark that some first tests using the scheme proposed in resulted in 
a similar performance gain as the two schemes presented above. 
In table |I| we give a typical example of the expense of a simulation following 
strategy (ii). We list both the cost of computing the lowest eigenvalue of Dj^D^ 
in terms of the number of iterations to minimize the Ritz functional |T9[ and the 
number of iterations to solve -D„X = t]. In both applications, a polynomial of 
degree n is used to approximate the square root. The numbers in table |I| indicate 
that a quenched calculation, employing Neuberger's operator, leads to a compu- 
tational cost that is comparable with a dynamical simulation using conventional 
operators. 



N 


n 




'^invert 


8 


190 


170 


80 


10 


250 


325 


200 


12 


325 


700 


300 



Table 1: is the number of lattice sites along a side of the hypercube; n, the 
degree of polynomial; n^vi the number of iterations required to obtain the lowest 
eigenvalue of Dj^Dn] and ninvert, the number of iterations necessary to compute 
X = 
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5 Conclusions 



The theoretical advance that an exact chiral symmetry brings to lattice gauge 
theory is accompanied by the substantial increase in numerical effort that is 
required to implement operators satisfying the GW relation. Thus, while the 
Nielsen-Ninomiya theorem has been circumvented, the "no free lunch theorem" 
has not. Whether alternative formulations, such as domain wall fermions, can 
help in this respect remains to be seen. 



Appendix 

We give here the explicit definitions needed in eq. (|T]). The forward and backward 
derivatives V^, V* act on a vector as 

V^$(x) = - [f/(x, /i)<l>(x + a/t) - $(x)] 

V*$(a;) = - [^{x) - U{x - ajl,fi)~^^{x - afi)] , 

where fi denotes the unit vector in direction /i. The (gauge) field U{x, /i) G SU (3) 
lives on the links connecting lattice points x and x + afi and acts on the colour 
index a = 1, 2, 3 of the field $. Finally, the Dirac matrices 7^, /i = 0, 1, 2, 3 are 
hermitean 4x4 matrices acting on the Dirac index i of the field $. Their explicit 
form is given by 



7m 





el 



(16) 



with 



Co = 1, Cfc = -iak 



(17) 



and 



01\ /0-i\ /lO 

"^='10 "^=U "^= 0-1 



(18) 
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With the choice of the 7 matrices given above, the matrix 75 = 70717273 is 
diagonal and given by 



We finally note that whenever repeated indices appear, they are summed over. 
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